Hazard in Firebaugh
Samantha Liu
2021-01-05
acs5_vars <- readRDS("acs5_vars.rds")
get_census_bg <- function(group) {
getCensus(
name = "acs/acs5",
vintage = 2019,
region = "block group:*",
regionin = "state:06+county:019",
vars = paste0("group(", group, ")")
) %>%
mutate(block_group = paste0(state, county, tract, block_group)) %>%
select(
!c(GEO_ID, state, county, tract, NAME) &
!ends_with(c("EA", "MA", "M"))
) %>%
pivot_longer(
ends_with("E"),
names_to = "variable",
values_to = "estimate"
) %>%
left_join(
acs5_vars %>%
select(name, label),
by = c("variable" = "name")
) %>%
select(-variable) %>%
filter(block_group == "060190084011")
}
get_census_tract <- function(group) {
getCensus(
name = "acs/acs5",
vintage = 2019,
region = "tract:*",
regionin = "state:06+county:019",
vars = paste0("group(", group, ")")
) %>%
mutate(tract = paste0(state, county, tract)) %>%
select(
!c(GEO_ID, state, county, NAME) &
!ends_with(c("EA", "MA", "M"))
) %>%
pivot_longer(
ends_with("E"),
names_to = "variable",
values_to = "estimate"
) %>%
left_join(
acs5_vars %>%
select(name, label),
by = c("variable" = "name")
) %>%
select(-variable) %>%
filter(tract == "06019008401")
}Building Data obtained from Geofabrik makes many components of OSM available for easy download. After clicking a few levels down, you can find a NorCal page with a 700MB file for download with these attributes. Another option is Microsoft’s AI-generated building footprints, which may be more or less accurate than OSM in different places. Their GitHub repo includes links to download the files as GeoJSONs by state, but California is massive, so a nice option is to use county-level disaggregations done by Patty Frontiera at Berkeley.
Hazard Data
Data Sources: Cal Adapt
Row numbers that result in a lexicon error: 46, 62, 108, 109, 134, 220
exheat_full <- readRDS("exheat_full.rds")
exheat_full %>%
group_by(year) %>%
summarize(
exheat_mean = mean(exheat, na.rm = T),
exheat_sd = sd(exheat, na.rm = T)
) %>%
ggplot(
aes(
x = year
)
) +
geom_ribbon(
aes(
ymin = exheat_mean - exheat_sd,
ymax = exheat_mean + exheat_sd
),
fill = "grey"
) +
geom_line(
aes(
y = exheat_mean
)
) +
labs(
x = "Year",
y = "Extreme heat days",
title = "Fresno County CBGs"
)
exheat_2020s <-
exheat_full %>%
filter(year >= "2020-01-01" & year <= "2029-12-31") %>%
group_by(cbg) %>%
summarize(exheat = mean(exheat, na.rm = T)) %>%
left_join(fresno_cbgs %>% dplyr::select(cbg = GEOID)) %>%
st_as_sf()
heat_pal <- colorNumeric(
palette = "Reds",
domain = exheat_2020s$exheat
)
leaflet() %>%
addProviderTiles(providers$CartoDB.Positron) %>%
addPolygons(
data = exheat_2020s,
fillColor = ~heat_pal(exheat),
fillOpacity = 0.5,
color = "white",
weight = 0.5,
label = ~paste0(cbg, ": ", exheat, " extreme heat days per year over the next decade"),
highlightOptions = highlightOptions(
weight = 2,
opacity = 1
)
) %>%
addLegend(
data = exheat_2020s,
pal = heat_pal,
values = ~exheat,
title = "Extreme heat days<br>per year over<br>next decade"
)Exposure Data
Equity Analys on Heat Exposure - Assumes extreme heat affects everyone in a given census tract the same way
bg_income <-
getCensus(
name = "acs/acs5",
vintage = 2019,
region = "block group:*",
regionin = "state:06+county:019",
vars = "group(B19001)"
) %>%
mutate(cbg = paste0(state, county, tract, block_group)) %>%
select(
!c(GEO_ID, state, county, tract, block_group, NAME) &
!ends_with(c("EA", "MA", "M"))
) %>%
pivot_longer(
ends_with("E"),
names_to = "variable",
values_to = "estimate"
) %>%
left_join(
acs5_vars %>%
select(name, label),
by = c("variable" = "name")
) %>%
select(-variable) %>%
separate(
label,
into = c(NA,NA,"income"),
sep = "!!"
) %>%
filter(!is.na(income)) %>%
mutate(
income = case_when(
income %in% c("Less than $10,000","$10,000 to $14,999","$15,000 to $19,999","$20,000 to $24,999") ~ "Less than $25,000",
income %in% c("$25,000 to $29,999","$30,000 to $34,999","$35,000 to $39,999","$40,000 to $44,999","$45,000 to $49,999") ~ "$25,000 to $49,999",
income %in% c("$50,000 to $59,999","$60,000 to $74,999") ~ "$50,000 to $74,999",
TRUE ~ income
)
) %>%
group_by(cbg, income) %>%
summarize(estimate = sum(estimate, na.rm = TRUE))Extreme Heat Day buckets 250-300 should be lumped with 225-250 because of how small the total number of households they affect
heat_exp_income <-
bg_income %>%
left_join(
exheat_full %>%
filter(year >= "2020-01-01" & year <= "2029-12-31") %>%
group_by(cbg) %>%
summarize(exheat = sum(exheat)),
by = c("cbg")
) %>%
filter(!is.na(exheat)) %>%
mutate(
exheat_tier =
case_when(
exheat < 175 ~ "150-175",
exheat < 200 ~ "175-200",
exheat < 225 ~ "200-225",
# exheat < 250 ~ "225-250",
# exheat < 275 ~ "250-275",
TRUE ~ "225-300"
)
) %>%
group_by(income, exheat_tier) %>%
summarize(estimate = sum(estimate, na.rm = T)) %>%
mutate(
income = income %>% factor(levels = rev(c("Less than $25,000", "$25,000 to $49,999", "$50,000 to $74,999", "$75,000 to $99,999", "$100,000 to $124,999", "$125,000 to $149,999", "$150,000 to $199,999", "$200,000 or more")))
)heat_exp_income %>%
rbind(
heat_exp_income %>%
group_by(income) %>%
summarize(estimate = sum(estimate)) %>%
mutate(exheat_tier = "Total")
) %>%
ggplot() +
geom_bar(
aes(
x = exheat_tier %>% factor(levels = c("150-175", "175-200", "200-225", "225-300", "Total")),
y = estimate,
fill = income
),
stat = "identity",
position = "stack"
) +
labs(
x = "Days of Extreme Heat",
y = "Number of households",
title = "Extreme Heat Day Exposure by Income",
subtitle = "Fresno County - cumulated over 10 years",
fill = "Household Income"
) +
coord_flip() +
theme(
legend.position = "bottom",
legend.direction = "vertical"
)
heat_exp_income %>%
rbind(
heat_exp_income %>%
group_by(income) %>%
summarize(estimate = sum(estimate)) %>%
mutate(exheat_tier = "Total")
) %>%
ggplot() +
geom_bar(
aes(
x = exheat_tier %>% factor(levels = c("150-175", "175-200", "200-225", "225-300", "Total")),
y = estimate,
fill = income
),
stat = "identity",
position = "fill"
) +
labs(
x = "Days of Extreme Heat",
y = "Number of households",
title = "Extreme Heat Day Exposure by Income",
subtitle = "Fresno County - cumulated over 10 years",
fill = "Household Income"
) +
coord_flip() +
theme(
legend.position = "bottom",
legend.direction = "vertical"
)
Difference-in-Difference
ACS Demographics
Total number of households recorded in Firebaugh according to ACS 2015-2019: 2321 Total number of households recorded in neighborhood of interest according to ACS 2015-2019: 844 Total number recorded given by FHA 2020: 2033
# housing costs as % of income
fb_housing_costs <-
get_census_tract("B25106") %>%
separate(
label,
into = c(NA, NA, "tenure", "income", "% of income"),
sep = "!!"
) %>%
filter(!is.na(`% of income`))
# HH Income - cbg level
fb_hhincome <-
get_census_bg("B19001") %>%
separate(
label,
into = c(NA, NA, "income"),
sep = "!!"
) %>%
filter(!is.na(income))
# aggregate vehicles available - cbg level
fb_vehicles <-
get_census_bg("B25046") %>%
separate(
label,
into = c(NA, NA, "tenure"),
sep = "!!"
) %>%
filter(!is.na(`tenure`)) %>%
mutate(
hh = c(342, 502), # obtained from fb_beds
veh_per_hh = estimate / hh
)
# bedrooms - cbg level
fb_beds <-
get_census_bg("B25042") %>%
separate(
label,
into = c(NA, NA, "tenure", "beds"),
sep = "!!"
) %>%
filter(!is.na(beds))
# units per structure - cbg level
fb_units <-
get_census_bg("B25032") %>%
separate(
label,
into = c(NA, NA, "tenure", "units"),
sep = "!!"
) %>%
filter(!is.na(units))
# hhsize by units in structure
fb_units_hhsize <-
get_census_tract("B25124") %>%
separate(
label,
into = c(NA, NA, "tenure", "HH Size", "units"),
sep = "!!"
) %>%
filter(!is.na(units))
# total pop by units in structure - cbg level
fb_units_pop <-
get_census_bg("B25033") %>%
separate(
label,
into = c(NA, NA, "tenure", "units"),
sep = "!!"
) %>%
filter(!is.na(units))
# heating fuel
fb_fuel <-
get_census_tract("B25117") %>%
separate(
label,
into = c(NA, NA, "tenure", "fuel"),
sep = "!!"
) %>%
filter(!is.na(fuel))
# computer and internet access by educational attainment - cbg level
fb_internet_education <-
get_census_bg("B28006") %>%
separate(
label,
into = c(NA, NA, "education", "computer", "internet"),
sep = "!!"
) %>%
mutate(
internet = ifelse(computer == "No computer", "No computer", internet)
) %>%
filter(!is.na(internet))
# internet access by income - cbg level
fb_internet_income <-
get_census_bg("B28004") %>%
separate(
label,
into = c(NA, NA, "income", "internet"),
sep = "!!"
) %>%
filter(!is.na(internet))
# median hhincome by hhsize
fb_hhsize <-
get_census_tract("B19019") %>%
separate(
label,
into = c(NA, NA, "HH Size"),
sep = "!!"
)
# hhtype by hhsize - cbg level
fb_hhtype_hhsize <-
get_census_bg("B11016") %>%
separate(
label,
into = c(NA, NA, "hh type", "hhsize"),
sep = "!!"
) %>%
filter(!is.na(`hh type`))
# householder age by income - cbg level
fb_age <-
get_census_bg("B19037") %>%
separate(
label,
into = c(NA, NA, "age", "income"),
sep = "!!"
) %>%
filter(!is.na(`age`))fb_housing_costs %>%
group_by(tract, tenure, income) %>%
summarize(estimate = sum(estimate, na.rm = TRUE)) %>%
mutate(
income =
income %>%
substr(1, nchar(.) - 1) %>%
factor(
levels =
c("Less than $20,000", "$20,000 to $34,999", "$35,000 to $49,999", "$50,000 to $74,999", "$75,000 or more")
)
) %>%
ggplot() +
geom_bar(
aes(
x = income,
y = estimate,
fill = tenure
),
stat = "identity",
position = "fill"
) +
labs(
x = "Household Income",
y = "Proportion of Households",
title = "Tenure by Household Income",
subtitle = "Firebaugh (Tract level)",
fill = "Tenure Classification"
) +
coord_flip() +
theme(
legend.position = "bottom",
legend.direction = "vertical"
)
fb_housing_costs %>%
group_by(tract, `% of income`, income) %>%
summarize(estimate = sum(estimate, na.rm = TRUE)) %>%
mutate(
income =
income %>%
substr(1, nchar(.) - 1) %>%
factor(
levels =
c("Less than $20,000", "$20,000 to $34,999", "$35,000 to $49,999", "$50,000 to $74,999", "$75,000 or more")
),
`% of income` = `% of income` %>% factor(levels = c("Less than 20 percent", "20 to 29 percent", "30 percent or more"))
) %>%
ggplot() +
geom_bar(
aes(
x = income,
y = estimate,
fill = `% of income`
),
stat = "identity",
position = "fill"
) +
labs(
x = "Household Income",
y = "Proportion of Households",
title = "Housing Costs as Percent of Income",
subtitle = "Firebaugh (Tract level)",
fill = "Percent of Income"
) +
coord_flip() +
theme(
legend.position = "bottom",
legend.direction = "vertical"
)
fb_housing_costs %>%
group_by(tract, `% of income`, income) %>%
summarize(estimate = sum(estimate, na.rm = TRUE)) %>%
mutate(
income =
income %>%
substr(1, nchar(.) - 1) %>%
factor(
levels =
c("Less than $20,000", "$20,000 to $34,999", "$35,000 to $49,999", "$50,000 to $74,999", "$75,000 or more")
),
`% of income` = `% of income` %>% factor(levels = c("Less than 20 percent", "20 to 29 percent", "30 percent or more"))
) %>%
ggplot() +
geom_bar(
aes(
x = income,
y = estimate,
fill = `% of income`
),
stat = "identity",
position = "stack"
) +
labs(
x = "Household Income",
y = "Proportion of Households",
title = "Housing Costs as Percent of Income",
subtitle = "Firebaugh (Tract level)",
fill = "Percent of Income"
) +
coord_flip() +
theme(
legend.position = "bottom",
legend.direction = "vertical"
)
fb_units_hhsize %>%
group_by(tract, tenure, `HH Size`) %>%
summarize(estimate = sum(estimate, na.rm = TRUE)) %>%
ggplot() +
geom_bar(
aes(
x = `HH Size`,
y = estimate,
fill = tenure
),
stat = "identity",
position = "fill"
) +
labs(
x = "Household Size",
y = "Proportion of Households",
title = "Tenure by Household Size",
subtitle = "Firebaugh (Tract level)",
fill = "Tenure Classification"
) +
coord_flip() +
theme(
legend.position = "bottom",
legend.direction = "vertical"
)
fb_units_hhsize %>%
group_by(tract, units, `HH Size`) %>%
summarize(estimate = sum(estimate, na.rm = TRUE)) %>%
mutate(units = ifelse(units == "1, detached or attached:", "1, detached or attached", units)) %>%
ggplot() +
geom_bar(
aes(
x = `HH Size`,
y = estimate,
fill = units
),
stat = "identity",
position = "fill"
) +
labs(
x = "Household Size",
y = "Proportion of Households",
title = "Units in Structure by Household Size",
subtitle = "Firebaugh (Tract level)",
fill = "Total Units in Structure"
) +
coord_flip() +
theme(
legend.position = "bottom",
legend.direction = "vertical"
)
fb_fuel %>%
group_by(tract, tenure, fuel) %>%
summarize(estimate = sum(estimate, na.rm = TRUE)) %>%
ggplot() +
geom_bar(
aes(
x = fuel,
y = estimate,
fill = tenure
),
stat = "identity",
position = "fill"
) +
labs(
x = "Fuel Source",
y = "Proportion of Households",
title = "Tenure by Fuel Source",
subtitle = "Firebaugh (Tract level)",
fill = "Tenure Classification"
) +
coord_flip() +
theme(
legend.position = "bottom",
legend.direction = "vertical"
)
fb_hhsize %>%
select(-tract) %>%
rename("HH Income" = "estimate") %>%
filter(!is.na(`HH Size`)) %>%
knitr::kable() %>%
kableExtra::kable_styling()| HH Income | HH Size |
|---|---|
| 23314 | 1-person households |
| 34167 | 2-person households |
| 31296 | 3-person households |
| -666666666 | 4-person households |
| 67533 | 5-person households |
| 44228 | 6-person households |
| -666666666 | 7-or-more-person households |
Median Income Bracket for the neighborhood is $30,000 - $34,999 Approximate number of vehicles per household for owned-homes: 2.62 Approximate number of vehicles per household for rented-homes: 0.92
fb_beds %>%
group_by(block_group, tenure, beds) %>%
summarize(estimate = sum(estimate, na.rm = TRUE)) %>%
ggplot() +
geom_bar(
aes(
x = beds,
y = estimate,
fill = tenure
),
stat = "identity",
position = "stack"
) +
labs(
x = "Number of Bedrooms",
y = "Proportion of Households",
title = "Tenure by Number of Bedrooms",
subtitle = "Specific neighborhood focus (CBG level)",
fill = "Tenure Classification"
) +
coord_flip() +
theme(
legend.position = "bottom",
legend.direction = "vertical"
)
fb_units %>%
group_by(block_group, tenure, units) %>%
summarize(estimate = sum(estimate, na.rm = TRUE)) %>%
ggplot() +
geom_bar(
aes(
x = units,
y = estimate,
fill = tenure
),
stat = "identity",
position = "stack"
) +
labs(
x = "Total Units in Structure",
y = "Total Households",
title = "Tenure by Total Units in Structure",
subtitle = "Specific neighborhood focus (CBG level)",
fill = "Tenure Classification"
) +
coord_flip() +
theme(
legend.position = "bottom",
legend.direction = "vertical"
)
fb_units_pop %>%
group_by(block_group, tenure, units) %>%
summarize(estimate = sum(estimate, na.rm = TRUE)) %>%
ggplot() +
geom_bar(
aes(
x = units,
y = estimate,
fill = tenure
),
stat = "identity",
position = "stack"
) +
labs(
x = "Total Units in Structure",
y = "Total Number of People",
title = "Tenure by Total Units in Structure",
subtitle = "Specific neighborhood focus (CBG level)",
fill = "Tenure Classification"
) +
coord_flip() +
theme(
legend.position = "bottom",
legend.direction = "vertical"
)
fb_internet_education %>%
mutate(
internet = ifelse(internet == "No computer", "", internet),
comp_int = paste(computer, internet)
) %>%
ggplot() +
geom_bar(
aes(
x = education,
y = estimate,
fill = comp_int
),
stat = "identity",
position = "fill"
) +
labs(
x = "Educational Attainment",
y = "Total Households",
title = "Computer and Internet Access by Educational Attainment",
subtitle = "Specific neighborhood focus (CBG level)",
fill = "Computer and Internet Access"
) +
coord_flip() +
theme(
legend.position = "bottom",
legend.direction = "vertical"
)
fb_internet_income %>%
mutate(
income =
income %>%
substr(1, nchar(.) - 1) %>%
factor(
levels =
c("Less than $10,000", "$10,000 to $19,999","$20,000 to $34,999", "$35,000 to $49,999", "$50,000 to $74,999", "$75,000 or more")
)
) %>%
ggplot() +
geom_bar(
aes(
x = income,
y = estimate,
fill = internet
),
stat = "identity",
position = "stack"
) +
labs(
x = "Household Income",
y = "Total Households",
title = "Internet Access by Household Income",
subtitle = "Specific neighborhood focus (CBG level)",
fill = "Internet Access"
) +
coord_flip() +
theme(
legend.position = "bottom",
legend.direction = "vertical"
)
fb_hhtype_hhsize %>%
filter(!is.na(hhsize)) %>%
group_by(block_group, `hh type`, hhsize) %>%
summarize(estimate = sum(estimate, na.rm = TRUE)) %>%
ggplot() +
geom_bar(
aes(
x = hhsize,
y = estimate,
fill = `hh type`
),
stat = "identity",
position = "stack"
) +
labs(
x = "Household Size",
y = "Total Households",
title = "Household Type by Household Size",
subtitle = "Specific neighborhood focus (CBG level)",
fill = "Household Type"
) +
coord_flip() +
theme(
legend.position = "bottom",
legend.direction = "vertical"
)
fb_age %>%
filter(!is.na(income)) %>%
mutate(
income = case_when(
income %in% c("$10,000 to $14,999","$15,000 to $19,999") ~ "$10,000 to $19,999",
income %in% c("$20,000 to $24,999", "$25,000 to $29,999", "$30,000 to $34,999") ~ "$20,000 to $34,999",
income %in% c("$35,000 to $39,999","$40,000 to $44,999","$45,000 to $49,999") ~ "$35,000 to $49,999",
income %in% c("$50,000 to $59,999","$60,000 to $74,999") ~ "$50,000 to $74,999",
income %in% c("$75,000 to $99,999", "$100,000 to $124,999", "$125,000 to $149,999", "$150,000 to $199,999", "$200,000 or more") ~ "$75,000 or more",
TRUE ~ income
)
) %>%
group_by(block_group, age, income) %>%
summarize(estimate = sum(estimate, na.rm = TRUE)) %>%
mutate(
income =
income %>%
factor(
levels =
c("Less than $10,000", "$10,000 to $19,999","$20,000 to $34,999", "$35,000 to $49,999", "$50,000 to $74,999", "$75,000 or more")
)
) %>%
ggplot() +
geom_bar(
aes(
x = income,
y = estimate,
fill = age
),
stat = "identity",
position = "fill"
) +
labs(
x = "Household Income",
y = "Proportion of Households",
title = "Householder Age by Household Income",
subtitle = "Specific neighborhood focus (CBG level)",
fill = "Age of Householder"
) +
coord_flip() +
theme(
legend.position = "bottom",
legend.direction = "vertical"
)
## Land-Use
leaflet() %>%
# addTiles(group = "OSM Base Map") %>%
addProviderTiles(providers$CartoDB.Positron) %>%
addPolygons(
data = fb_parcels,
color = "green",
fillOpacity = 0.5,
weight = 1,
label = ~APN
) %>%
addPolygons(
data = fb_osm_bldg,
color = "#e851a2",
fillOpacity = 0.5,
weight = 1,
label = ~ID
)projection <- "+proj=utm +zone=10 +ellps=GRS80 +datum=NAD83 +units=ft +no_defs"
bldg_osm <-
fb_osm_bldg %>%
st_transform(projection)
parcels <-
fb_parcels %>%
st_transform(projection)
outside_area <-
bldg_osm %>%
mutate(bldg_area = st_area(.)) %>%
arrange(desc(bldg_area)) %>%
st_difference(st_union(parcels))
yards <-
parcels %>%
mutate(lot_area = st_area(.)) %>%
arrange(desc(lot_area)) %>%
st_difference(st_union(bldg_osm))
leaflet() %>%
addTiles() %>%
addPolygons(
data = yards %>% st_transform(4326),
fillColor = "blue",
color = "white",
weight = 1,
fillOpacity = 0.5
)blocks <-
parcels %>%
st_union()
street_edges <-
blocks %>%
st_difference(blocks %>% st_buffer(-1)) %>%
st_intersection(parcels, .) %>%
filter(APN %in% yards$APN)
leaflet() %>%
addTiles() %>%
addPolygons(
data = yards %>% st_transform(4326),
fillColor = "yellow",
color = "white",
weight = 1,
fillOpacity = 0.5
) %>%
addPolygons(
data = street_edges %>% st_transform(4326),
fillColor = "brown",
color = "brown",
weight = 2
)bldg_join <-
bldg_osm %>%
st_centroid() %>%
st_join(parcels) %>%
st_set_geometry(NULL) %>%
left_join(bldg_osm %>% select(ID)) %>%
st_as_sf()
street_edges_buffered <-
1:nrow(street_edges) %>%
map_dfr(function(x){
temp_apn <-
street_edges[x,] %>%
pull(APN)
temp_bldg <-
bldg_join %>%
filter(APN == temp_apn) %>%
st_union()
if (length(temp_bldg) == 0) { return(street_edges[x, ]) }
nearest_distance <-
st_nearest_points(street_edges[x, ], temp_bldg) %>%
st_length()
street_edges[x, ] %>%
st_buffer(nearest_distance)
})
leaflet() %>%
addTiles() %>%
addPolygons(
data = yards %>%
st_transform(4326),
fillColor = "yellow",
color = "white",
weight = 1,
fillOpacity = 0.5
) %>%
addPolygons(
data = street_edges_buffered %>%
st_transform(4326),
fillColor = "brown",
color = "brown",
weight = 2
)yards_setbacks <-
parcels %>%
mutate(lot_area = st_area(.)) %>%
arrange(desc(lot_area)) %>%
st_buffer(-4) %>%
st_difference(st_union(bldg_osm %>% st_buffer(6))) %>%
st_difference(st_union(street_edges_buffered))
leaflet() %>%
addTiles() %>%
addPolygons(
data = yards_setbacks %>% st_transform(4326),
fillColor = "blue",
color = "white",
weight = 1,
fillOpacity = 0.5
)buildable_area <-
yards_setbacks %>%
st_buffer(-3.9) %>%
st_buffer(3.9, joinStyle = "BEVEL") %>%
filter(!st_is_empty(.)) %>%
st_cast("MULTIPOLYGON") %>%
st_cast("POLYGON") %>%
mutate(
buildable_area = st_area(.) %>% as.numeric()
) %>%
filter(buildable_area >= 160)
leaflet() %>%
addTiles() %>%
addPolygons(
data = yards_setbacks %>%
st_transform(4326),
fillColor = "red",
color = "white",
weight = 1,
fillOpacity = 0.5
) %>%
addPolygons(
data = buildable_area %>%
st_transform(4326),
fillColor = "green",
color = "white",
weight = 1,
fillOpacity = 1
)Avg Buildable Area per Single-Family DU:
avg <- ((buildable_area %>% st_union() %>% st_area()) - (outside_area %>% st_union() %>% st_area())) / nrow(fb_osm_bldg)
avg## 935.0226 [ft^2]Scale-up to the entire city of Firebaugh:
fb_sfdu <-
get_census_tract("B25032") %>%
separate(
label,
into = c(NA, NA, "tenure", "units"),
sep = "!!"
) %>%
filter(units == "1, detached")
avg * sum(fb_sfdu$estimate)## 1297811 [ft^2]Linear Programming
max_units <- function(units, area, people, b1, b2, b3, b4) {
# Set coefficients of the objective function
f.obj <- c(1, 1, 1, 1)
# Set matrix corresponding to coefficients of constraints by rows
# Do not consider the non-negative constraint; it is automatically assumed
f.con <- matrix(c(1, 1, 1, 1,
700, 900, 1200, 1500,
2, 4, 6, 8,
1, 0, 0, 0,
0, 1, 0, 0,
0, 0, 1, 0,
0, 0, 0, 1
), nrow = 7, byrow = TRUE)
# Set inequality signs
f.dir <- c("<=",
"<=",
"<=",
"==",
">=",
">=",
">="
)
# Set right hand side coefficients
f.rhs <- c(units,
area,
people,
b1,
0.31 * (b2 + b3),
0.43 * (b2 + b3),
b4
)
# Final value (z)
print(lp("max", f.obj, f.con, f.dir, f.rhs))
# Variables final values
lp("max", f.obj, f.con, f.dir, f.rhs)$solution
}result <- max_units(844, 345080.1, 1469, 73, 19, 294, 14)## Success: the objective function is 322.455sum(result * c(2, 4, 6, 8))## [1] 1469Electric Usage
electric <-
map_dfr(paste0("Q", 1:4), function(x) {
if (x == "Q4") { year <- 2019 } else { year <- 2020 }
file <-
read_csv(paste0("PGE_", year, "_", x, "_ElectricUsageByZip.csv")) %>%
filter(ZIPCODE == 93622) %>%
filter(CUSTOMERCLASS == "Elec- Residential")
})sum(electric$AVERAGEKWH) * 0.31 ## [1] 2196.97Gas Usage
gas <-
map_dfr(paste0("Q", 1:4), function(x) {
if (x == "Q4") { year <- 2019 } else { year <- 2020 }
file <-
read_csv(paste0("PGE_", year, "_", x, "_GasUsageByZip.csv")) %>%
filter(ZIPCODE == 93622) %>%
filter(CUSTOMERCLASS == "Gas- Residential")
})sum(gas$AVERAGETHM) * 52.85## [1] 18867.45